These exercises require the latest version of popdemo. Head to popdemo’s GitHub (github.com/iainmstott/popdemo) for installation instructions (don’t forget also to load the package using library(popdemo)!)

Complete the core exercises (in normal print) first, as code in each section continues from the last. Afterward, return to the “extras” sections (in italics), and try writing your own code. Code to be run is in chunks:

# a comment
an_input()
##   an output

Key terms in the text are in bold italic, and functions or arguments are fixed width.

1. Data

We will use a matrix projection model (MPM) to explore population dynamics for the desert tortoise Gopherus agassizzii, with medium fecundity1. The population is found in the Mojave desert, USA. There are 8 stages are based on age and size (carapace length in mm):
- Yearling (age 0-1)
- Juvenile 1 (<60 mm)
- Juvenile 2 (90-99mm)
- Immature 1 (100-139mm)
- Immature 2 (140-179mm)
- Subadult (180-207mm)
- Adult 1 (208-239mm)
- Adult 2 (>240mm)

Load in the data:

The numbers in the matrix (called matrix elements or transitions) describe the probability of moving FROM stages in each column TO stages in each row, within the time interval chosen. For example, for this desert tortoise matrix, in any year a subadult (stage 6) has approimately 24.9% probability of becoming an adult (stage 7): this may be called a growth or progression transition. Likewise, in any year a subadult has about 67.8% chance of staying a subadult: this is called a stasis transition. This means that 100 - (67.8 + 24.9) = 7.3% of subadults die every year. There are different types of transitions: in this matrix there are also fecundity transitions which describe offspring production, and subadults produce on average 1.3 offspring per year. Other species may have different transitions, including skipping stages through fast growth, shrinkage or fission (especially in modular organisms, e.g. most plants, corals), or asexual reproduction.

Matrix elements combine underlying vital rates such as survival, growth and reproduction. For example, the subadult to adult transition (24.9%) combines probability of growing to adult size in one year, and probability of surviving the year. The subadult fecundity (1.3) combines the average number of offspring produced by subadults per year, and probability that those offspring survive the year.
 

A baby desert tortoise

2. Deterministic projections

The core function for understanding population dynamics is the project function. It can be used to understand a number of different types of population dynamics, including deterministic models, stochastic models, long-term dynamics and short-term dynamics.

In a deterministic model, there is no density dependence, and no stochasticity: the vital rates of the model, and therefore the matrix elements, don’t change from timestep to timestep (in a density dependent model, they would change according to the size of all or part of the population, and in a stochastic model they would change randomly at each iteration of the model to incorporate random environmental or demographic variation). The projected population dynamics (time series of population size and structure over time) come from multiplying the matrix and the starting poulation vector:

\[\textbf{n}_\text{t} = \textbf{A}^t\textbf{n}_0\]

That is, the population size at time t nt is equal to projection matrix A multiplied by the population size at time 0, n0. We have the matrix, but not a vector. We will:
- choose a vector using a random uniform distribution
- project this vector using the project function from popdemo
- plot the projection

The project function should usually have a matrix passed to the A argument, and a vector, which can be a single vector, multiple column vectors in a matrix, "n" (a default option which projects a set of stage-biased vectors) or "diri" (which projects large set of random vectors). We will encounter all these options in this vignette. The time argument gives the number of projection intervals, but the output will have length time + 1, because the population size at time 0 is included.

The project function returns an object of class ‘Projection’, containing the overall population size over time. In this case, this is a single vector but for objects containing multiple projections (as we’ll encounter later), the projections are held in a matrix with one column per projection (therefore the number of rows is equal to time + 1).

The ‘Projection’ object also includes information on the time series of population vectors (i.e. (st)age-specific number / density). For an object containing a single projection, each row represents one timestep and each column represents one (st)age. We will access vectors for time intervals from 0 to 10:

The time series of stage 2 sizes is:

When a Projection object contains multiple projections, the vectors are contained in a 3-dimensional array. As for a single projection, the first dimension is time and the second is (st)age, whilst the third represents each separate population projection. Therefore, for example, the first 10 timesteps of the third projection would be accessed using object[1:11, ,3]; the sizes of the second stage for all projections would be accessed using object[ , 2, ]; the size of the second stage at the 10th timestep of the second stage in the third projection would be object[1:11, 2, 3].

Population growth in an MPM is geometric: when you plot population size on a log scale (as we have here), it’s easy to see that the population settles to a stable geometric rate of decline. At this point, the population also has a fixed structure: the relative numbers of individuals in each stage don’t change (although the absolute numbers do; this is what causes the population decline). These stable long-term dynamics are often called asymptotic dynamics. The short-term dynamics that happen before this stable state is reached and are different to asymptotic dynamics are called transient dynamics. In this exercise we’ll explore both asymptotic and transient dynamics of deterministic MPM models.
 

ii. EXTRAS…
’Projection objects contain further information that can be accessed in similar ways to the vec() function. Type ?Projection-class to see further options.

Try altering the parameters in the project function: change the amount of time the model is projected using the time argument, or change the population vector using the vector argument. Have a look at how the asymptotic dynamics don’t change as the vector changes, but observe how the transient dynamics change. See whether changing the vector changes the amount of time taken to reach stable state. Type ?project to see further options (some of them will be covered in this vignette).

The plot function takes any of the usual graphical parameters: try changing the lines to points (type argument), changing the line colour, type or thickness(col, lty, lwd), changing the box around the plot (bty), or any other graphical parameters (see ?par). There are further options available specifically for ‘Projection’ objects: see ?Projection-plots.
 

3. Asymptotic dynamics

In the long term, the population dynamics are described by the dominant eigendata of the matrix:

3.1 Asymptotic growth

The eigs function returns the dominant eigendata (or “eigenstuff”) of the matrix2. The growth rate is commonly referred to as lambda (\(\lambda\)). In this case, \(\lambda\) < 1 which means the population declines. If \(\lambda\) = 1 the population neither grows nor declines, and If \(\lambda\) > 1, the population grows.

3.2 Asymptotic population structure

The rest of the eigenstuff describes other aspects of stable dynamics. “ss” refers to the stable structure: this is the ratio of numbers of individuals in each stage once the population reaches stable state, and the vector is usually denoted with w. A population with this structure, then it will grow/decline at the stable rate from the outset (populations starting with a stable structure will always be shown with dashed lines):

w is scaled so that ||w||1 = 1.3

3.3 Reproductive value

“rv” refers to the reproductive value vector, and is usually denoted with v. This is the contribution that each stage makes to stable growth (through survival, growth and reproduction). Stages with high current and reproduction and survival have high reproductive value. v is scaled so that vTw = 14 when ||w||1 = 1.

 

iii. EXTRAS…
The eigs function allows you to choose what eigenstuff you want to calculate. Try replacing “all” with one or more of “lambda”, “ss”, or “rv”. If you want to look at subdominant eigenstuff then the base function “eigen” does this, with eigen(A) giving the right eigenvectors and eigen(t(A)) giving the left eigenvectors.
 

4. Transient dynamics

Before settling to stable growth rate, a population will grow, decline or fluctuate in growth at rates faster and/or slower than asymptotic growth. We can see this in our plot of population dynamics of the desert tortoise. These population dynamics are called transient dynamics and are a little harder to characterise than asymptotic dynamics as they’re so variable and depend on the population structure5.

4.1 Standardisations

popdemo contains functions that calculate deviations from stable growth at various points along the population projection. These transient indices make two standardisations: starting population vector n0 is scaled by ||n0||1 so that it sums to 1:

\[\hat{\textbf{n}}_0 = \frac{\textbf{n}_0}{||\textbf{n}_0||_1}\]

The projection matrix is scaled by \(\lambda\), so that lambda of the scaled matrix becomes 1:

\[\hat{\textbf{A}} = \frac{\textbf{A}}{\lambda}\]

These standardisations allow comparison of transient dynamics between populations with different sizes, and with different long-term dynamics. We can visualise this this in a standardised population projection:

Transient dynamics vary according to the starting population vector. If there is an overrepresentation of individuals in stages with high survival and/or fertility, then the population will grow faster (or decline slower) than the stable rate (this is called ‘amplification’). If there is an overrepresentation of individuals in stages with low survival and zero/low fertility, then the population will grow slower (or decline faster) than the stable rate (this is called ‘attenuation’). Indices of amplification are always >1, which means that log-transformed indices of amplification are always >0. Indices of attenuation are always <1 but >0. This means that log-transformed indices of attenuation are always <0.

4.2 Transient indices

We can measure transient density at any time along the projection (we say “density” because a standardised dynamic is no longer directly equivalent to size… but this is not the same thing as spatial density!). But in comparative analysis, to make things comparable between populations, we can use comparable timepoints.

Two possibilities are at t = 1 and at \(t\to\infty\). These indices are called reactivity and inertia respectively:

4.2.1 Reactivity

Reactivity (r1) is the population size in the first timestep of the projection (one year), relative to a population with stable growth. A reactivity of 2 would mean that in the first timestep, the population grows twice as fast as its stable growth rate. Reactivity is calculated using:

\[\textit{reactivity} = P_1 = ||\hat{\textbf{A}} \hat{\textbf{n}}_0||_1\]

As we can see, measurements of transient dynamics use the standardisations we encountered earlier. These “standardisations” mean that transient indices are measured relative to long-term population growth, and initial population size. The interpretation of this is that transient indices measure the ratio of population size compared to a population growing at a stable rate.

4.2.2 Inertia

Inertia (i1) is the ratio of the size of the population in the long-term, relative to a population with stable growth. An inertia of 4 would mean that after the transient period, the population settles to a size 4 times as large as a population that grows with stable rate. Inertia is calculated using:

\[\textit{inertia} = P_\infty = \frac{\textbf{v}^\text{T}\hat{\textbf{n}}_0||\textbf{w}||_1}{\textbf{v}^\text{T}\textbf{w}}\]

w and v are the right and left eigenvectors as described in the “Asymptotic dynamics” section above. If w and v are scaled as described, then \(\textit{inertia} = \textbf{v}^\text{T}\hat{\textbf{n}}_0\).

4.2.3 Maximum amplification and attenuation

Comparable points at which to measure transient dynamics don’t have to be at the same time of the projection. maximum amplification and maximum attenuation describe the largest and smallest population sizes (relative to long-term growth \(\lambda\)). Maximum amplification and maximum attenuation can occur at any point along the projection. Use return.t = TRUE to return the time at which they occur, as well as their value.

Some populations only have a maximum amplification (because they never attenuate), some only have a maximum attenuation (because they never amplify), and some have both. With this proviso in mind, these indices are calculated using: \[maximum \ amplification = \bar{P}_{max} = \max_{t > 0}||\hat{\textbf{A}}^t\hat{\textbf{n}}_0||_1\] \[maximum \ attenuation = \underline{P}_{min} = \min_{t > 0}||\hat{\textbf{A}}^t\hat{\textbf{n}}_0||_1\]

If we want to compare several different population structures, then it’s possible to project several simultaneously. We will:
- create two extra vectors; one which amplifies and one which attenuates
- bind these with the random vector into a matrix with 3 columns
- project the model with this matrix passed to the vec argument
- calculate reactivity, inertia and maximum amplification / attenuation for the extra vectors
- plot these on the graph

 

iv. EXTRAS…
Some further functionality offered by the reac, inertia, maxamp and maxatt functions. Try using, for example, return.N = TRUE to give the transient population size (including influences of initial population size and asymptotic population growth), and use this to plot transient indices on non-standardised population projections.
 

5. Transient bounds

Transient dynamics are very variable, and there are an infinite number of different starting population vectors. Sometimes we don’t know what the population structure is (making a fair census of plants and animals is difficult!), but it is possible to get an idea of what the transient dynamics will be like. Transient bounds capture the outer extremes of transient dynamics; all population trajectories lie within the bounds.

5.1 Plotting transient bounds

It’s easy to plot the bounds by adding bounds = TRUE when plotting a population projection:

5.1.1 Stage-biased projections

The bounds are calculated from the set of projections of “stage-biased” vectors, each of which has 100% of individuals in one stage. For example, the stage-biased vectors for a 3-by-3 matrix are [1,0,0], [0,1,0] and [0,0,1]. To project these, use vector = "n"; this is the default value, so in fact if vector isn’t set then the stage-biased vectors will be projected automatically:

The top and bottom lines are the bounds on population dynamics. No deterministic projection can be outside these lines.

5.2 Bounds on transient indices

When calculating transient indices, use bound = "upper" or bound = "lower" to calculate bounds on reactivity and inertia. To calculate bounds on maximum amplification or attenuation, just don’t pass anything to vector:

Using \(\rho_1\) to refer to reactivity bounds, \(\rho_\infty\) to refer to inertia bounds, \(\rho_{max}\) to refer to the maxmimum amplification bound, \(\rho_{min}\) to refer to the maximum attenuation bound, and with overbars to refer to upper bounds and underbars to refer to lower bounds6:

\[\bar{\rho}_1 = ||\hat{\textbf{A}}||_1 \hspace{0.75cm} and \hspace{0.75cm} \underline{\rho}_1 = minCS(\hat{\textbf{A}})\] \[\bar{\rho}_\infty = \frac{v_{max}||\textbf{w}||_1}{\textbf{v}^\text{T}\textbf{w}} \hspace{0.75cm} and \hspace{0.75cm} \underline{\rho}_\infty = \frac{v_{min}||\textbf{w}||_1}{\textbf{v}^\text{T}\textbf{w}}.\] \[\bar{\rho}_{max} = \max_{t > 0}||\hat{\textbf{A}}^t||_1 \hspace{0.75cm} and \hspace{0.75cm} \underline{\rho}_{min} = \min_{t > 0}(minCS(\hat{\textbf{A}}^t))\]

When w and v are scaled as decribed in the “Asymptotic dynamics” section, then \(\bar{\rho}_\infty = v_{max}\) and \(\bar{\rho}_\infty = v_{min}\).
 

v. EXTRAS…
We’ve seen a little bit of the diversity of different transient dynamics that can be shown by a model. But there are an infinite number of different possible starting population vectors. If population vectors are drawn at random, it is possible to shade in the space within the transient bounds to work out the likelihood of different transient densities. This can be done by using vector = "diri" in the project function, which draws populations at random from a dirichlet distribution and projects them. Plot this using plottype = "shady":

Dirichlet projections have their own options within the project and plot functions: see the “Projection-class” and “project” help pages.

  1. Doak et al. (1994) Ecol. Appl., 4, 446-460.↩

  2. for a detailed introduction to matrix models and eigendata, see Caswell (2001) Matrix Population Models, Sinauer.↩

  3. ||w||1 is the one-norm of w, which is equal to its sum.↩

  4. vT is the transpose of v: not to be confused with an exponent!↩

  5. Stott et al. (2011) Ecol. Lett., 14, 959-970 contains a review on measuring transient dynamics in MPMs.↩

  6. \(||\hat{\textbf{A}}||_1\) is the one-norm of \(\hat{\textbf{A}}\), is equal to its maximum column sum. \(minCS\) describes the minimum column sum. \(v_{max}\) and \(v_{min}\) are the maximum and minimum entries of v, respectively.↩